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Two-level system with a thermally fluctuating transfer matrix element: Application to 

the problem of DNA charge transfer 

Maria R. D'Orsogna and Joseph Rudnick 
Physics Department, University of California, Los Angeles, CA 90095-1547 
(Dated: February 1, 2008) 

Charge transfer along the base-pair stack in DNA is modeled in terms of thermally-assisted 
tunneling between adjacent base pairs. Central to our approach is the notion that tunneling between 
fluctuating pairs is rate limited by the requirement of their optimal alignment. We focus on this 
aspect of the process by modeling two adjacent base pairs in terms of a classical damped oscillator 
subject to thermal fluctuations as described by a Fokker-Planck equation. We find that the process 
is characterized by two time scales, a result that is in accord with experimental findings. 

PACS numbers: 87.15.-v, 73.50.-h,82.30.Fi 



I. INTRODUCTION 

In spite of the fact that a decade has passed since the 
first definitive observation of charge transfer along the 
DNA base-pair stack jjj , the detailed properties of this 
process have not been definitively elucidated. This is 
partly due to the inherent complexity of the molecular 
structure of DNA, and to the large number of exter- 
nal and intrinsic factors that exert an influence on DNA 
structure and behavior. The current unsettled situation 
also reflects the absence of an overall agreement on the 
precise mechanism by which this charge transport takes 
place. One of the key issues that awaits full illumina- 
tion is the role of disorder — both static and dynamic — 
on the propagation of charge along the base-pair stack. 
A related, and quite fundamental, question is whether 
charge transport is a coherent quantum mechanical pro- 
cess, like conduction of electronic charge against a static, 
or deformable, background, or whether it it takes place as 
fundamentally incoherent transport, as a variation of the 
random walk. The answers to these and other questions 
will have a significant impact on both our understand- 
ing of the biological impact of charge transport in DNA 
and the development of applications based on this phe- 
nomenon. 

Despite the often contradictory results of experimen- 
tal investigations j|, ||, ||, ||, a few conclusions seem 
inescapable. The first is that long-range charge trans- 
port along the base-pair stack depends quite strongly on 
the sequence of the base pairs Q. In addition, base-pair 
mismatches can have a significant deleterious effect on 
charge transport ||, |^| (see, however ]l0|]). Furthermore, 
strands of DNA display considerable disorder, both static 
and dynamic |§ [u| p). Finally, several sets 
of experiments on ensembles of short DNA strands have 
uncovered an unusual two-step charge transfer process 
p6| , [l7| . These studies focus on fluorescent charge donors 
intercalated in DNA oligostrands. As the charge migrates 
towards the acceptor, the fluorescence is quenched and 
the rate of migration is determined by the decaying flu- 
orescence profile. The data reveals that this decay pro- 
cess occurs according to two characteristic time-scales 



which are separated by more than an order of magni- 
tude jl6|, [l7| . Any model that purports to explain charge 
transport must take all this into account. 

In this paper, we discuss a model for short range 
charge transport along a base pair stack that undergoes 
substantial structural fluctuations. The process occurs 
via thermally-assisted quantum mechanical tunneling of 
charge carriers from one base pair to the next, under 
the assumption that this tunneling is properly character- 
ized as occurring in the presence of a dissipative environ- 
ment. A key conjecture is that charge transfer takes place 
only when the neighboring pairs are in a state of optimal 
"alignment" , and that this alignment is statistically un- 
likely in thermodynamic equilibrium. As we will see, this 
conjecture leads in a natural way to a model exhibiting 
the dual-time-scale feature described above. Addition- 
ally, the model generates predictions that can be readily 
tested. We shall relate the problem at hand to the dy- 
namics of a simple two level system (TLS), realized by a 
donor and an acceptor state. 

In Section |J, we briefly recapitulate what is known 
about the tunneling process in the presence of friction for 
a TLS system. We also quantify our notion of a coordi- 
nate 9 associated with the "alignment" of adjacent base 
pairs and of the influence of the dynamics of this new 
coordinate on charge transfer. Section III specifies the 
model for describing a generic collection of two-level sys- 
tems (TLS), initially in the donor state and characterized 
by a fluctuating alignment variable 0. The probability 
distribution of donor states, W(9, 9, t) obeys a Kramers 
equation with a sink term due to charge transfer to the 
acceptor. The rate of charge transfer will be expressed 
by the fluctuating rate r(#). This Kramers expression is 
recast into the form of a Volterra equation with the use 
of a Lie- Algebra approach defined on the Hilbert space 
of the eigenstates of the Kramers equation for T(9) = 0. 
We will discuss limiting cases of the solution to obtain 
physical insight and to reveal the two-time-scale decay 
of the probability distribution due to the sink term. We 
conclude in Section |iy| with a discussion of the possible 
application of our results to charge transfer in strands of 
DNA consisting of several base pairs. 

The key result of our calculations lies in the determi- 
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nation of P(9*,t), the probability distribution of donor 
states evaluated at the optimal configuration 9* and with 
the 9 variable integrated out. Indeed, under the assump- 
tion that the tunneling process is most effective at 9 ~ 9* , 
this quantity is directly related to the fluorescence in- 
tensity I(t) of the base pair complexes, as probed by J. 
Barton and A. Zewail [|l6| |l7), through the following: 



I(t) = Jo 



i-r / p{o*,t')dt' 

Jo 



(1) 



The quantity Iq of the above relationship is a propor- 
tionality constant and T is the integrated rate of transfer 
to the acceptor. We shall determine the double expo- 
nential character of P(9*,t), and hence of I(t), in qual- 
itative agreement with the experimental findings. The 
conjectures made on the existence of an optimal and un- 
likely configuration 9* will be crucial in obtaining the 
two stage decay process, a result that justifies the as- 
sumptions made. 

The model we shall construct is obviously not re- 
stricted in application to DNA oligostrands. Using our 
results, we may conclude that in an ensemble of generic 
systems the migration of a particle from donor to ac- 
ceptor proceeds statistically as a two-time scale process, 
provided the transfering process is of rare occurrence. 



II. THE TUNNELING PROCESS 

The process of charge transfer from a donor site to an 
acceptor site — a two-level system — is ubiquitous in bio- 
chemical and physical phenomena |l8|. It occurs under 
a broad variety of spatio-temporal conditions. Chemical 
bond formation or destruction, ATP production in pho- 
tosynthetic reactions, or the operation of semiconduct- 
ing devices, all involve the transfering of charges to and 
from specific sites, via thermal activation or quantum- 
mechanical tunneling through an energy barrier. Be- 
cause of its intrinsic nature, charge transfer via quantum- 
mechanical tunneling takes place on a length scale of up 
to tens of angstroms 19 1; larger distances are possible 
if other transport mechanisms are involved. These in- 
clude thermal hopping among sites, which are typical in 
disordered systems, the creation of conduction bands in 
metals, or of lattice distortions of polaronic type in spe- 
cific systems. 

Quantum-mechanical tunneling from a donor site to 
an acceptor site is quite simply represented by a two- 
level system (TLS) [gOj. In this description, the tunneling 
particle is limited to being in the donor or in the acceptor 
state, while the other degrees of freedom of the system, 
nuclear for instance, describe the charge potential energy. 

The energetic profile of the system is thus characterized 
by a multidimensional surface of which the acceptor and 
the donor states constitute relative minima, separated by 
a barrier. Of the many existing degrees of freedom, it is 
often possible to identify a "reaction coordinate" y such 



V(y) 




FIG. 1: This figure illustrates the nature of the tunneling 
transition. The two parabolic curves shown correspond to 
the two versions of the potential V(y,a z ) in Eq. (0), one 
corresponding to the "donor" state in which the tunneling 
particle is on one site and the other to the "acceptor" state, 
in which the particle is on the other one. The two energies 
Ef and Et, = Ef — e referred to in the text are the forward 
and backwards barrier energies respectively. The horizontal 
axis corresponds to the reaction coordinate, y. 



that the energy barrier between donor and acceptor is 
minimized along this specific direction. The progress of 
the reaction is then dominated by the evolution along 
this coordinate and the potential energy surface can be 
reduced to an effective one-dimensional curve. 

In certain systems the physical interpretation of the 
reaction coordinate is immediate: it may be the relative 
bond length in two diatomic molecules, or solvent polar- 
ization around the donors and acceptors ]2l] j. It is not 
an easy task to give a physical interpretation of the re- 
action coordinate in the case of DNA base pairs because 
of the many possibilities involved - intra-base distance, 
mobile counter-ion concentration, solvent concentration, 
or a combination of all the above. A possibility is offered 



by Ref. 22 where it is suggested that the most relevant 
quantity is the interaction of the charge with the polar 
water molecules of the solvent. In this paper we shall 
refer to the reaction coordinate y in most general terms. 

A common representation of tunneling with dissipation 
is through the spin-boson formalism |20| . The donor and 
acceptor states are represented by means of a pseudo 
spin, which points up when the charge is in the donor 
state and down otherwise. The Hamiltonian of the sys- 
tem is given by: 



where 



H E T = Ta x + -J- + V(y,a z )+H ha , th , (2) 



V(y,<J z ) = ^MLo 2 (y + yo<J z ) 2 + \ea z , (3) 
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and <j XtZ are the Pauli matrices. The charge in the donor 
(up) state corresponds to the potential V(y, +) whose 
equilibrium reaction coordinate is — j/o, aud the converse 
state corresponds to V(y,—), whose stable minimum is 
at yo- The i?bath term represents contributions to the 
Hamiltonian of a dissipative environment coupled to the 
reaction coordinate. Figure [l] illustrates the meaning 
of the potential V(y,a z ) in the effective Hamiltonian of 
Eq.(||). The curve marked A corresponds to the poten- 
tial term in the donor state, while the curve marked B 
represents the potential function in the acceptor state. 

This model has been thoroughly analyzed in the work 
by Garg et al. S based on earlier work by Leggett pi[ . 
A similar analysis, but within a more chemical frame- 
work, is presented by Marcus et al. |25| ]. Energy con- 
servation requires that charge transfer takes place only 
when the reaction coordinate is close to the degeneracy 
point y = y* for which V(y*,+) = V(y*,— ); once the 
degeneracy point is reached, charge transfer is possible 
only because of the non zero off-diagonal tunneling ma- 
trix elements r. 

The tunneling rate T from donor to acceptor, is calcu- 
lated in the above references. For moderate dissipation 
of the reaction coordinate, it is given by: 



h \E r k B T cS 



e —Ef/kBT B tt _|_ e ~E b /kBT B {{ 



(4) 



where the reorganization energy E r and the energy bar- 
riers Ef and Eb depend on the details of the potential de- 
scribed by the reaction coordinate. In the limit of high 
temperatures T e g reduces to the usual temperature T, 
whereas in the opposite limit the quantity is tempera- 
ture independent. 

The novelty explored in this paper is the introduction 
and investigation of the effect of a second reaction co- 
ordinate, 9, governing the charge transfer process and 
coupled not to the energy, but to the off-diagonal tunnel- 
ing element r, hitherto been treated as a constant, and 
which we now write t(8). 

This new coordinate reflects the conjecture that in the 
case of DNA the tunneling matrix element is highly sensi- 
tive to the donor-acceptor relative configuration. Charge 
transport along DNA in fact occurs along the stacked 
base pairs by means of overlapping n orbitals, and at 
room temperature, these base pairs strongly fluctuate 
with respect to each other through variations of the twist, 
tilt and roll parameters [2(J . The existence of base pair 
fluctuations for DNA in solution is very well established, 
and is corroborated by experimental |p_ 2|| and molecular 

For such a highly asym- 



dynamics studies [^3[ [l4| 
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metric system such as DNA, fluctuations in the relative 
orientation of donors and acceptors affect the magnitude 
of the orbital overlap between pairs, and the new collec- 
tive coordinate 9 embodies the effects of these fluctua- 
tions. 



We will also assume that the 9 variable is slowly vary- 
ing compared to the motion of the reaction coordinate y, 
so as to define the lowest energy scale of the system. We 
may then separate the motion of the two reaction coordi- 
nates in a Born-Oppcnhcimcr spirit. Charge transfer will 
be assumed to be instantaneous once the optimal 6 = 6* 
value is reached, and a purely classical framework will be 
utilized for the ^-dynamics. The new reaction coordinate 
9 need not necessarily be pictured as a geometrical one, 
although this is the framework we will be utilizing in this 
paper. Just as in the case of the y reaction coordinate, 9 
may be associated to the particular chemical environment 
of the molecule or to any other quantity influencing the 
strength of the tunneling element t between the donor 
and the acceptor sites. 

The Hamiltonian describing the system thus, is a mod- 
ified version of the spin-boson Hamiltonian introduced in 
Eq. (||) with a t(9)<j x off-diagonal term, as also described 
in earlier work |27f| . In order for charge transfer to take 
place, we will assume that the reaction coordinate cou- 
pled to the energy must be close to the degeneracy point 
y = y* , and, also, that the 9 coordinate must be in the 
neighborhood of an optimal value 9*, which maximizes 
the tunneling amplitude. The physical picture to asso- 
ciate to this requirement is that the relative "alignment" , 
9, does not favor charge transfer unless an optimal con- 
figuration is reached: t(9) ~ unless 9 ~ 6* . This 
conjecture will prove to be crucial in yielding the two 
time-scale charge transfer of references [|l6| [l7| . 

In analogy to the experimental work cited above, we 
consider a collection of such two- level systems, with the 
charge initially located on the donor site. Each one 
of these systems is associated to a particular t(9) and 
through Eq. (||) to a particular T(9) rate. Our objective 
is to determine the mechanisms of charge transfer taking 
into account the 9 time evolution and the T(6) rates ac- 
cordingly distributed. We shall assume the 9 dynamics 
to be governed by small, Langevin type random fluctua- 
tions. At t = 0, when the external charge is injected on 
the donor site, the distribution of 9 values is the usual 
Boltzmann distribution. If the occurrence of the optimal 
9* configuration is relatively unlikely, we will indeed be 
able to show that the transfer process is characterized by 
a two time scale migration of the initial donor population. 

The emergence of two time scales in the transfer pro- 
cess can be physically explained as follows. The exis- 
tence of an initial non-zero population of TLS presenting 
the optimal value 9* , ensures that rapid tunneling to the 
acceptor. The 9 distribution is thus depleted of popu- 
lation at the special value and other transitions are for- 
bidden to take place. The other TLS will tunnel to the 
acceptor only after the system has re-equilibrated and re- 
populated the optimal configuration, a process which is 
slow, because of the assumption that the optimal config- 
uration is a relatively unlikely one. Hence, the existence 
of a fast, initial decay followed by a slower decay process. 
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III. THE TLS AND 9 FLUCTUATIONS 



ables: 



The model 



t') = 2q5(t-t'). (6) 



Consider a collection of TLS which at the initial time 
t = are all in the up-donor configuration, and charac- 
terized by the angular parameter 9. Let us denote by 
W(9, 9, t) the TLS population remaining in the up-donor 
state at time t and for which the collective angular vari- 
able and its velocity are specified. 

The physical requirement that 9 be randomly, classi- 
cally, fluctuating in time, translates into the fact that 
W(9, 9, t) must evolve according to a Fokker-Planck type 
equation as dictated by standard Langevin theory. To 
this probability evolution equation we must add an addi- 
tional depleting term, that which represents tunneling to 
the donor site as given by the T(9) term discussed above. 

Different scenarios are possible for the 9 dependence of 
t and hence of T. As discussed in the above section we 
shall focus on the particular situation in which tunneling 
is possible only for a very specific subset of energetically 
unfavorable 9 values. In this picture, tunneling is al- 
lowed only if donors and acceptors reach an optimal — but 
unlikely — orientation one with respect to the other. By 
including the tunneling term in the time evolution equa- 
tion for W(9, 9, t) we obtain a modified Fokker-Planck 
equation that may be used to approach any physical sys- 
tem in which the presence of a depleting term competes 
with the usual Langevin fluctuations. The most natural 
choice for the 9 motion, the one we shall discuss in the 
remainder of this paper, is that of a damped harmonic 
oscillator. We shall see that starting from an initially 
equilibrated system in which the 9 distribution is the the 
Boltzmann one, the insertion of the tunneling term will 
result in the emergence of the two time scales discussed 
above. We will refer to the time derivative of the 9 coor- 
dinate as u. The rotational moment of inertia associated 
to 9 is denoted by / and its rotational frequency by tt. 

The goal of the next subsections will be to determine 
W(9* , u, t), and in particular its integration with respect 
to the u variable. As described in the introduction in 
fact, it is this quantity that is directly related to the 
experiments we wish to model by means of Eq. dll) . 



B. Kramers equation with a sink term 

The generic damped harmonic oscillator subject to ran- 
dom noise responds to the following Langevin-type equa- 
tions: 



- 1 u-tt 2 9 + i 1 {t), 



(5) 



where the stochastic force r](t) is assumed to be a zero- 
mean gaussian and whose correlation function is dictated 
by the fluctuation-dissipation theorem for classical vari- 



The corresponding Fokker-Planck equation may be writ- 
ten by identifying |2^] the proper coefficients in the 
Kramers-Moyal expansion from Eq. (|^) and is generally 
referred to as the Kramers equation. This equation gov- 
erns the time evolution of the distribution, W(9,u,t), of 
an ensemble of systems obeying the equations of motion 
(||). It takes the form: 



dW 



dW d_ 
l ~d9 + !hi 



[{ lu + n 2 9)w\ 



d 2 W 
' du 2 



(7) 



The above equation is thoroughly analyzed in |30| ] , where 
assuming an initial probability distribution W(9,u,0) = 
5(9 — 9')5(u — u'), the probability W(9,u,t), as well as 
other relevant statistical quantities, are obtained. At 
equilibrium Kramers equation is solved by the time inde- 
pendent Boltzmann distribution, W(9, u, t) — ipo.o(9,u) 
with: 



7$1 

ipo,o(9,u) = - — exp 
2nq 



7 I 2 
2q V 



fi 2 



(8) 



Under the assumptions discussed earlier, the probability 
distribution function W (9, u, t) for a particle localized on 
the donor site and describing an effective angle 9 with its 
neighbor, will be described by the time evolution equa- 
tion for a collection of damped oscillators subject to a 
decay term T, representing tunneling to the acceptor. 
The latter term is appreciable only for a specific value of 
the 9 coordinate 9*: 



dW 
~dt 



= HW-T(6,u,t) W. 



(9) 



The H term is the differential operator that stems from 
the right hand side of Eq. (0). We shall assume the 
decay term to be introduced at time t = 0, prior to which 
the system had attained its equilibration state. In other 
words, we choose the initial distribution W(9, u, 0) to be 
Boltzmann-like, as expressed in Eq. (||). For simplicity, 
we choose T(9) to be independent of u and of t and to be 
a gaussian centered on 9* and with width a: 



T(9) 



exp 



l*\2 



2a 



(10) 



The coefficient k contains the physical parameters of tem- 
perature and energy as expressed in Eq. (Q). We also 
impose the constraint that at t = the optimal value 
9* carries a small Boltzmann weight. This is equivalent 
to the physical assumption that the occurrence of par- 
ticle tunneling is a rather unlikely event, and that the 
system tends to relax to 9 values that are far from the 
tunneling point. We also impose the width of the decay 
gaussian ^/a, to be small compared to 9*, so that T(9) 
is highly pea ked aro und the optimal configuration value 
9*: % /7o : < y/q~/W^ 9*- 
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In the following subsections we will solve Eq. (§) for 
the early and long time regimes. The general solution 
for arbitrary times is contained in the appendix. The 
coupling of the system to the orientational degree of free- 
dom along the lines discussed above, manifests itself very 
clearly in the unusual time dependence of the probability 
distribution. Two different decay rates in fact arise, with 
a rapid initial decay of the donor population W{9, u, t) 
followed by a slower transfer process. The ratio of these 
two time scales, and the main result of this analysis is 
succinctly expressed by Eq. ( |27| ) in terms of all the phys- 
ical parameters of this system. 



Such a solution will eventually exceed the "leading or- 
der" one. We determine the coefficients A m ^ n by re- 
quiring that there be no secular term in the solution to 
the equation. It is precisely this constraint that con- 
stitutes the underlying idea of multiple scale analysis. 
The above condition translates into requiring that the 
non-homogeneous term within parenthesis in Eq. (|l^) or 
equivalently in Eq. (fla) vanish: 



(16) 



C. Short time regime 

In order to determine the asymptotic behavior of 
W(6,u,t) in the early time regime, we consider Eq. (^) 
with the gaussian choice of T(9) and we perform a multi- 
ple time scale analysis ]3l| . This is carried out by intro- 
ducing a new ad-hoc variable £ = T(9)t, into the proba- 
bility distribution, and by seeking solutions in the form 
W(9,u,t) = W (6,u,t,t)+T(9) W 1 (9,u,t,£) + .... The 
Fokker-Planck equation is thus expanded in powers of 
T(9) and, for the zeroth and first order terms, it yields: 



We now solve for A m ^ n . Imposing the initial condition 
W(0,u, 0) = ^,0(9,11) and reinserting £ = T(9)t the so- 
lution reads: 



W (9, u, t) = Vo,o(0, u) exp [-T{9)t]. (17) 

The above is a zero-th order approximation to the full 
problem presented in (12) to the extent that the effect 
of H acting on tT{9) can be neglected with respect to 
T(9) itself. In other words, Eq. Jl7j ) is an approximate 
solution as long as: 



dW 
dt 



HW = 0, 



HWi 



dW 



Wo 



(11) 

uT w Wl{l2) 



Note that the partial derivative with respect to t in the 
above equations treats £ as an independent variable. The 
solution to the first equation is expanded in terms of the 
complete set of functions ^ m .n(9, i t, t) that solve Eq. jll| ) 
- obtained in Eq. (Al) and Eq. (A£) of the appendix - 
with coefficients A m , n that depend on £, i.e: 



(9, u) e 



A m „t 



(13) 



t < 



r(0) 

\uT e (9)\ \u(9~9*)\' 



(18) 



This equation is valid only under the conditions expressed 
in dl8| ) and up to t ~ T~ 1 (9). For this time limitation to 
be meaningful, it is necessary that the width of the decay 
term ^fa be finite. In the limit that the width vanishes 
the above analysis fails, since the expansion parameter 
diverges. At time t ~ we cannot approximate T{9) 
by a strict delta function. Note that for 9^9*, the 
tunneling point, and for finite u the condition arising 
from the multiple scale analysis (t ~ \j2-na / ' n) is the 
most stringent one, and the probability distribution is 
approximated by: 



Substituting this solution for Wq into Eq. (Ji2|), the in- 
homogeneous term in square brackets becomes: 



W(6*,u,t) =^M0*,lt) exp [- 



Kt 



(19) 



(14) 



If this were the only term present on the right hand side 
of Eq. (|l2|), then Wi(9, u, t, £) would contain a secular 
term in its solution of the type: 



We now perform an i nteg ration over the u variable 
on both sides of Eq. (|17|) and obtain an approxima- 
tion for the distribution probability function P(9, t) = 
J^W&u,*) du: 



P(6,t)~M6) exp [-r(0) t]. 



(20) 



Wi(6,u,t) 



(15) 



where V'o(^) is the Boltzmann distribution associated to 
the 9 variable ipo{9) = ipofi{6,u) du. For small 
times, P(9,t) retains its initial gaussian shape, with its 
amplitude decreasing exponentially. 
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D. Long time regime 

In this subsection we determine the long time asymp- 
totic behavior of W(9,u,t), utilizing some of the results 
obtained in the appendix for arbitrary tim es. In particu- 



l ar, w e adapt the kernel expansion of Eq. ( A15 ) an d Eq . 
( [A 16] ) to the long time regime. Differentiating Eq. ( |A15| ) 
with respect to t and with the gaussian choice for T(9) 
we obtain: 

^ = -JJ dO'du' ^ {9' ,u')T{9') 



K(9,9',u,u',0) W(9',u',t) + 

' dt ' W^' 9 '' ^ w ^' 7 u '' 1 ~ ^ 



(21) 



where the integrals in 9' and in u range from — oo to +oo. 
The time-derivative of the kernel in the last integral can 
be obtained with the use of the expression obtained in 



Eq. ( A16) but with the summation restricted to non-zero 



values of the integers m and n. The contribution to the 
kernel of the term associated with m = n = is time- 
independent, and it has the form i/}Q t o(9,u) i^o t o{9' , u'). 
We then replace dt> K with dt>K' where K' is defined as 
the kernel without the first (m, n = 0) summand. 

The function K' and its time derivative contain expo- 
nentially vanishing terms in t. The time integrand in Eq. 
( pl| ) will therefore be appreciable only for t' < ft- 1 where 
ft c is a cutoff frequency of the order of | Ax,o j = ft. For t » 
ft- 1 we can approximate W(9',u',t — t') ~ W{9',u',t) 
and restrict the time interval from the origin to ft" 1 . In- 
tegrating by parts, and using the above approximation 
for W(9', u',i), the time integral yields: 



dW 



dO'du' V^o (#',0 r(6»') 
iw{6',u',t) K(9,9',u,u',0)+ 

'(9, 9'u, u', ft" 1 ) - K'(9, 9', u, u', 0) } 



(22) 



K 



This equality is simplified by K'(9, 9'u, u', ft^ 1 ) being 
negligible. We can now rewrite the right hand side of 
Eq. (|22] ) as: 



dW 



d9' du' 



%l{9',u')T{9>) (23) 



M9,u) Tp 0fi (e',u') W(9',u',t) 



Since we are dealing with non-zero times, the 9' integra- 
tion can be performed under the assumption that T(9') 
is highly peaked around 9* and T(9) ~ k 5(9 — 9*): 

-^ = -«Vo,oM J du' W(9*,u',t). (24) 

A last integration in the u variable, performed on both 
sides of the equation, yields the probability distribution 



function for the 9 variable 

dP(9,t) 
dt ~ ~ 



kM8) P(9*,t). 



(25) 



For 9 — 9* the above relationship yields a decay rate of 
— ft ipo(9*), and for arbitrary 9 values we obtain the its 
behavior in the late time regime: 

P(0,t) = P MO) exp [-ft Vo(0*) t] . (26) 



E. The two time scales 

As anticipated, two different scenarios for P(9* , t) emerge 
from the a naly sis carried out in the previous subsections. 
From Eq. (|20|), at early times, the decay to the acceptor 
state is rapid, occurring at a rate ri = ft/\/27r<7, whereas 
at latter times the rate is as given above: r-i — nipo(9*). 
The ratio between the two is 



r-2 



k B T 
alfl 2 



exp 



/ft 2 
2fc B T 



[9*f 



» 1, 



(27) 



as follows from the assumptions made on the gaussian 
T(9). The initial decay is much faster than that at later 
times. 



Numerical results 



Based on the general solution of Eq. ( A15[ ), we 
present a numerical analysis of the distribution function 
W(9,u,t) for different choices of its arguments. In this 
equation the probability distribution W(9, u, t) is cast in 
a Volterra-type formulation, for which solutions can be 
constructed iteratively in time. The probability distribu- 
tion W(9, u, t) as expressed in Eq. ( A15) in fact, depends 



only on its previous history and on the known propagator 
function. 

For a numerical approach, it is necessary to dis- 
cretize the 9,u,t variables and keep track of the value of 
W(9,u,t) for every position and velocity at every tem- 
poral iteration. While feasible, this approach is rather 
cumbersome, since for every time step — kAt we must 
create a new 0(N 2 ) matrix W(9i,uj,tk), 1 < i,j < N, 
where N is the number of spacings for the position and 
velocity meshes. On the other hand, the evaluation at of 
W(9*, Uj, tk) where 9* represents the 9i interval centered 
on the optimal value 9* is greatly simplified if the corre- 
sponding mesh is chosen so that T(9) may be replaced for 
all purposes by a delta function at non-zero times. The 
recursive equations now involve only the O(N) element 
vector W(9*, Uj , t k ), 1 < j < N. 

At t = 0, when the propagator itself is a point source, 
the gaussian shape for T(9) must be retained for finitc- 
ness, but the iteration at a time that is far from zero does 
not involve values of the position that are significantly 
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different from 9* . The it mesh is chosen with Ait = 0.05 
and the time interval spacing is At = 0.01. 
In order to insure consistency with the constraint yfd <C 
\J qj^l^l 2 < 0* we choose the following parameters: a = 
10~ 4 , jft 2 = 2q, 9* = 1.5. The a parameter for the 
underdamped case is chosen as a = 0.02, whereas k is 
fixed at k — 0.4. The resulting probability distribution 
W(9* ,u,t) is plotted in Figure [| as a function of u for 
various time intervals. 

Two features of the evolving distribution are notewor- 
thy. The first is the depression around u = 0. The second 
is a clear asymmetry in the velocity distribution, in that 
the distribution for negative values of the velocity, u, is 
lower than for positive u values. The reason for the first 
feature is the fact that when the velocity is low, a pair 
will remain in a nearly optimal configuration longer, and 
hence a tunneling event, leading to depletion of the dis- 
tribution, is more likely. The asymmetry can be ascribed 
to the fact that the optimal orientation is at positive 
values of the parameter 9. The time evolution equation 
encapsulates two mechanisms, one pushing the distribu- 
tion towards its Boltzmann limit, the other being the 
tunneling process that leads to depletion of the distribu- 
tion at values of 9 close to 9* . In light of the trajectory of 
the underdamped oscillation, a member of the ensemble 
with negative velocity, u, is likely to be within a half an 
oscillation period of having passed with a small velocity 
through 9* , which is positive, while a representative with 
positive u is more likely to have spent more then half an 
oscillation period away from the optimal tunneling con- 
figuration. This latter, positive u configuration will have 
had more time to experience the "restorative" effects of 
the mechanism that acts to generate the Boltzmann dis- 
tribution. 

It is also possible to perform a it-variable integration 
and obtain the time dependence of P(9*,t). The param- 
eters are chosen as above, and the two time scale decay 
of P(9*,t) can be clearly seen to occur with rates n and 
T2 as described in Eq. (|27]) . Also note that both at large 
and short times P(9*,t) is proportional to W(9*,0,t). 
The above results, and the expressions for r\ and are 
not affected by changes in the damping variable a. As 
anticipated, Figure [J clearly shows the double exponen- 
tial decay of P(9* , f ), in agreement with the experimental 
results of ||, 0. 

IV. DISCUSSION 

The model we have presented is expected to be of sig- 
nificant relevance to charge transfer in DNA. Thermal 
fluctuations strongly affect the structure of molecule, and 
an accurate description requires this motion to be taken 
into account. 

Not only has the existence of fluctuations been experi- 
mentally documented jl2| , but it has also been suggested 
flEf that the motion that most affects the electronic cou- 
pling between base pairs - what we have referred to as 



0.04 




-4.0 -2.0 0.0 2.0 4.0 
u (arb. units) 

FIG. 2: The probability distribution at the optimal configu- 
ration W(9*,u, t) for various time intervals. The top curve is 
evaluated at t = and is the initial Boltzmann distribution, 
evaluated at the unlikely configuration 9*. The remaining 
curves are its time evolution, up to t — 5 of the lower curve. 
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FIG. 3: P{6* ,t) for the parameters chosen in the text (dotted 
curves). Fhe inset pertains to early times. The solid curves 
are drawn for comparison and are exponential decays e~ rt , 
with a rate ri = /t/V27rcr for the short times of the inset, and 
r : 2 = Kipo{6*) for the long time regime. Note the two distinct 
time scales. 



t(9) - is their sliding one with respect to the other. It 
must be pointed out that both these studies focus on 
DNA in solution, not on dry strands of DNA. 

On the other hand, charge transport with more than 
one rate has been reported in the literature [fT6| . For 
an oligomer with the ethedium molecule acting as the 
donor, charge transfer is found to occur along the same 
patterns as described by our model, with two time scales 
of 5 and 75 picoseconds. Two-time-scale decays are also 
observed in a series of measurements |fl7| performed on 
shorter strands of donor and acceptor complexes (Ap-G). 
In these experiments the Ap donor can be treated, for all 
practical purposes, as an intrinsic purine base, and the 
ambiguity related to the choice of an extraneous donor 
(the ethedium of the previous reference) is removed. 

In both these experiments, an increase in the length 
results in a competition between the fast and slow ex- 
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ponential decays in favor of the slower time component. 
Increasing the length of the system diminishes the possi- 
bility that multiple base pairs simultaneously arrange in 
the configuration that facilitates rapid charge transfer. 
When the process of optimal alignment does occur (a 
relatively likely event only for a few base pairs), the tun- 
neling might not even require localization of the charge 
on each base pair, and super-exchange can take place. 

For long strands of DNA, thus, we expect the two in- 
trinisic rates associated to a single charge transfer to be 
averaged out in favor of the slower component. Traces of 
this unusual two time scale migration mechanism how- 
ever, may be found in the fact that DNA conductivity 
is enhanced upon increasing the temperature |32|], pre- 
sumably allowing for greater base pair motion. Charge 
transfer is also hindered by disruptions to the stacking, 
which alter the base pair's ability to find optimal trans- 
fer configurations, such as the insertion of bulges along 
the helix or of strong mismatches within the base pair 
stacking 33, [54| which are poorly compatible with the 
intrinsic conformation of the aromatic pairs. Lastly, it 
is noted that charge transfer effectiveness seems to be 
inversely proportional to the measured hypochromicity 
p5[ , a quantity that determines the ordering of base pairs 
along a certain direction and defined as the reduction of 
absorption intensity due to interactions between neigh- 
boring electric dipoles. From this data it is apparent that 
the higher the disorder of the system, the more efficient 
charge transfer is. It would be interesting to see how dif- 
ferent solvent environments affect conduction along the 
molecule in relation to their effect on structural fluctu- 
ations. More temperature-dependent experimental mea- 
sures are desirable as well. 



V. CONCLUSIONS 

We have presented a model for a spin boson TLS whose 
tunneling matrix element depends on the structural con- 
formation of the donor with respect to the acceptor. In 
the limit that the relative geometry between the two fluc- 
tuates in time defining the lowest energy scale, we are 
led to a classical problem, that of a collection of damped 
harmonic oscillators obeying a modified Fokker-Planck 
equation. If charge transfer proceeds only for specific 
orientations of the donor with respect to the acceptors, 
the resulting rate for charge transfer is divided into a fast 
component at short times and a subsequent slower one. 
These results agree with the experimental findings of two- 
time-scale charge transfer in the donor intercalated DNA 
complexes of J.Barton and coworkers [l7|. It must be 
noted that an implicit assumption of this work is that for 
long range DNA conduction mediated by thermal fluctu- 
ations once the charge has undergone a transfer between 
base pairs it does not return to the pair at which it is 
originally localized. However, it is reasonable to assume 
that the transfer process will continue after this event 
has occurred and that subsequent events will, with some 



probability, deposit the charge at its point of origin at 
a later time. We have performed calculations on a two- 
time-scale hopping model based on the results obtained 
here j36|. In these calculations, the single set of two 
base pairs is replaced by a linear array. We have de- 
termined the probability that the charge carrier is at its 
point of origin as a function of time, t, after its having 
been placed there. We find that this probability exhibits 
two-time-scale behavior, with an initial, brief, rapid, ex- 
ponential decay followed by a much slower, power-law, 
decay at later times times. The long-time asymptotics of 
this process are those of a random walk. 
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APPENDIX A: GENERAL SOLUTION OF THE 
KRAMERS EQUATION 

We shall adopt a Lie-Algebra approach |37| to identify 
a complete set of orthonormal functions that solve the 
homogeneous problem in the general case of Eq. ^ , and 
through them the general solution for the decay equation 
(0) will be found. 

Let us look for solutions of the following type, where 
m and n represent non negative integers: 



* m ,„(6>, u,t) = tp m ^ n (9,u)e 



(Al) 



Upon insertion of the above expression in Eq. (^) a time 
independent Schrodinger-like equation can be written: 



■(A m , n + i)i>m,n = H'lpr, 



where: 



H'(e,u) = qpl + 'Yup u + Q? 



upe, 



(A2) 



(A3) 



and the subscripts represent derivatives, p u = d/d u . As 
expected, the time independent Boltzmann distribution 
satisfies the homogeneous equation, as can be verified by 
direct substitution with Ao.o = 0. The physical require- 
ment that solutions must be well behaved as t — > oo, i.e. 
that the A m) „'s be non negative, suggest that this is the 
ground state: 



groun 



d(0,u,t) = tjjofl(9,u). 



(A4) 



The other solutions are found by constructing the ladder 
operators. For the underdamped case, we introduce the 
a variable such that cos a = 7/ (2k;) and impose that 
[H', O] = I O with I and O respectively complex variable 
and operator to be determined. In practice, the operator 
O corresponds to either a raising or a lowering operator. 
Two sets of solutions exist for the following 'quanta' ii j: 



li = fie 



(A5) 
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for which the associated raising and lowering operators 
Ri 2 and L\ 2 are: 



Pi, 2 = ~Pe + h,2 Pu, 

12 



(A6) 



Li, 2 = n l e + l Pg + i h2 (l Pu + u). (A7) 

The commutation rules for the above operators can be 
easily derived as: 

[Ri,Rj] = 0, [L^Lj] = 0, [R U L 2 ] = 0, 
[R 1 ,L 1 ] = n 2 (e- 2 * a (A8) 
[R 2 ,L 2 ] = n 2 (e+ 2la -l). 

The raising operators applied to the ground state yield 
the set of solutions ip m . n for Eq. (|a|) with the associated 
eigenvalues X m ,n as follows: 

Wn (0, «) = i?2 ^ V>0,0 (0, U), (A9) 

A TO: „ = TOfie" la + nfie iQ . (A10) 

It is worth noting that the Hamiltonian H' can also be 
reformulated as H' = (2Qi sin a)" 1 [(L 2 R 2 ) - (Z-iPi)]. 
In order to construct solutions to the non-homogeneous 
problem within the Hilbert space spanned by the set of 
solutions {ip m . n (0, u)}, it is necessary to determine the 
orthonormality of those solutions. To this purpose, let us 
consider the following W„ l:n (9,u)} = {P 2 n P^il) Ofi (0, «)} 
where Pi 2 are operators defined as: 



Pi, 



Zl.2 Pu- 



(All) 



We can now prove an orthogonal relation between the two 
sets, using the commutation rules and and introducing 
iPqq(6,u) as a weighting function: 



du dQ (f)' m , n ,(0,u) ip l(9,u) ^ TO) „(0,w) = 

^m,n 6 mrn ' 6 nn i. (A12) 



The integration limits are over the entire real axis, both 
for 9 and u. The orthonormal set of eigenfunctions is 
thus expressed as {C~^„ 4>' mn (0,u)}, to which we refer 
as {4> m ,n{0, u)}. The constant of proportionality C m ^ n is: 

^-j (l- e - 2 ») m (l-e 2M )rA13) 

Let us now look for the full solution W(9, u, t) to Eq. (||), 
posing it in the following form: 



W(6,u,t) =y^ j h m<n (t) ip m<n (6,u) e 



A m „t 



(A14) 



The h m , n (t) functions are to be determined, in analogy to 
the scattering problem of particles in quantum mechan- 
ics. Let us assume that the decay term is introduced 
at time t = 0, and that the initial distribution is the 



equilibrium solution to the homogeneous problem, i.e. 
t he g round state as expressed in Eq. (||). Inserting Eq. 
([All ) in Eq. (§) and using the orthonormality relations, 
it is possible to find time evolution equations for h m , n (t) 
and to write a recursion formula for the full solution: 



W(0,u,t) = W(0,u,O)- / dt' 

Jo 



<W / du' 

j J — OO 

K(9,9' 1 u,u' 1 t-t') ^l(0\v!) 



T(0',u',t') W(0',u',t') 



(A15) 



Here, we have kept T a generic function of all variables 
and the K function is the response kernel of the system: 



K(6,9',u,u',t) = 



(A16) 



The product W'(0,u,t) = K(0,9',u,u',t) ^(9' ,u'), is 
the distribution function for the homogeneous system, 
under the initial conditions W'(0, u, 0) = 6(0 — 0')6(u — 
it'). Its asymptotic behavior reduces to the Boltzmann 
distribution, and apart from t = 0, it is an analytical 
function in all its variables. The explicit representation 
of the kernel may be written by insert ing t he expressions 
for ip m ,n(0,u) and 4> m ,n(0,u) in Eq. (kit): 



K(0,u,0'u',t) 



(A17) 



exp 



exp 



g (d e - fle+ ia d u ) (d e> + ile +ia d u .) _ ne+l 
7^ 2 (1 - e+ 2la ) € 

q (d e - Qe-* a d u ) (dg> + Qe-* a d ul ) _ Qe - t 



ip ,o(9,u) Vo,o(0V), 

where the exponential terms are intended as operators 
acting on the ground state wave functions. As it is writ- 
ten, the above kernel is still expressed symbolically. In 
order to obtain its explicit form it will suffice to perform 



a Fourier transform of Eq. ( A17), and then return to real 



space, a straightforward but tedious calculation we omit. 
The complete solution for the kernel is given by [B8| : 



K(0,0',u,u',t) = (^ 



1 



(A18) 



exp 



exp 



-L-(n 2 (l-n)(9 + 9') 2 + (l + l)(u-u') 
+ 2mQ(0 + 0')(u-u' 



4qG 



( fl 2 (l + n)(0 - 9') 2 + (1 - l)(u + u') 
+ 2mQ(0' -0)(u + u') 



10 



In order to keep a lighter notation, we have suppressed 
the time dependence of the T(t), G(t), l(t), m(t), n(t) 
functions. They are defined as: 



l(t) sin a 
m{t) sin a 
n(t) sin a 



-ate 



,-ntc 



-ntc 



iQ sin(a + ntsina), (A19) 



sin(Ot sin a), 



(A20) 



sin(a — fit sin a). (A21) 



The functions T(t) and G(t) are combinations of the 
above: 



In order to ensure integrability for Eq. ( A15 ), some limi- 
tations are posed on the form of the T(0', u', t') function. 
For instance, the seemingly most natural choice, a delta 
function centered around 9* , yields a non integrable ex- 
pression for W{9, u, t) at small times, when the kernel 
is a product of delta functions itself. Instead, the gaus- 
sian choice introduced earlier, with its finite a, ensures 
integrability at all time regimes. 



T(t) = l + l{t) -n(t) -n(t)l{t)-m 2 {t), (A22) 
G{t) = l + n{t)-l{t)-n(i)l{t)-m 2 {t). (A23) 
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